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1.  INTRODUCTION  AND  PROBLEM  STATEMENT 

Micro  air  vehicles  (MAVs)  -  vehicles  with  a  maximum  dimension  of  no  more  than  15  cm  -  offer 
attractive  solutions  for  applications  such  as  surveillance,  reconnaissance,  and  search  and  rescue 
operations.  The  reduction  in  size  is  accompanied  by  a  change  in  aerodynamic  flight  regime, 
characterized  by  a  relatively  low  Reynolds  number.  In  this  regime,  flapping  flight  introduces  several 
aerodynamic  advantages,  as  evidenced  by  the  ubiquitous  use  of  flapping  flight  by  various  animals  (for 
example,  about  75%  of  the  approximately  13,000  species  of  warm-blooded  vertebrates  fly). 

Many  aspects  of  flapping  flight  aerodynamics,  however,  are  yet  to  be  fully  understood,  such  as 
unsteady  aerodynamic  lift-enhancing  phenomena,  identification  and  characterization  of  flow  separation 
and  reattachment  physics,  laminar-to-turbulent  transition,  and  transient  vortex  dynamics.  Animal  flight 
in  nature  is  replete  with  a  wide  variety  of  features,  the  physics  and  benefits  of  which  are  also  not  yet 
fully  understood,  such  as  compliant  and  adjustable  wings,  widely  varying  flapping  wing  kinematics,  and 
extreme  maneuverability. 

The  purpose  of  the  research  performed  under  this  grant  was  to  explore  aerodynamics  and  flow-solid 
interactions  of  a  flapping  wing  mechanism.  This  was  accomplished  by  performing  research  in  the 
following  four  areas: 

1.  Development  and  testing  of  a  flapping  mechanism  featuring  adjustable  wing  kinematics. 

2.  Exploration  of  different  kinematic  representations  of  flapping  wing  trajectories. 

3.  Study  of  hardware-in-the-loop  optimization  approaches. 

4.  Three-dimensional  flow  imaging  using  synthetic  aperture  particle  image  velocimetry  (SAPIV). 

The  outcomes  of  these  four  activities  are  summarized  below. 

2.  DEVELOPMENT  AND  TESTING  OF  FLAPPING  MECHANISM  WITH  ADJUSTABLE  KINEMATICS 

A  robotic  flapping  wing  mechanism  for  use  in  general  studies  involving  flapping  flight  and 
laboratory-based  experimental  optimization  of  flapping  trajectories  was  designed  and  fabricated  (see 
Fig.  1).  The  design  allows  for  dynamic  adjustment  of  flapping  trajectories  in  air  or  liquids  with  three 
rotational  degrees  of  freedom  on  each  wing.  The  mechanism  design  has  adjustable  kinematics,  can 
operate  in  water,  oil,  or  air,  is  fixed,  and  can  be  used  for  flow  visualization,  particle  image  velocimetry 
(PIV),  and  force  analysis.  The  mechanism  is  described  in  detail  in  [1,2].  The  mechanism  was  used  in  a 
custom-built  tank  (see  Fig.  2  and  Ref.  [3]). 
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Figure  1.  Renderings  of  flapping  mechanism. 
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Figure  2.  Flow  tank  schematics.  Left:  Cross 


-section  showing  submerged  mechanism.  Right:  Tank  only. 


3.  THEORETICAL  COMPARISON  OF  WING  TRAJECTORY  KINEMATICS 

In  the  work  described  in  [1,2],  the  parameters  of  a  simple  Fourier  series  representation  of  the  wing 
angles  were  optimized  using  a  hardware-in-the-loop  approach.  The  follow-on  study  summarized  here 
investigated  the  use  of  other  kinematic  representations.  The  objective  was  to  explore  the  relative 
strengths  of  the  candidate  representations,  refine  the  optimization  methods,  and  implement  the  most 
promising  in  hardware-in-the-loop  optimization  to  find  the  best  parameters.  A  novel  kinematic 
representation  was  developed  to  enable  modulation  of  the  down-stroke-to-up-stroke  ratio  of  a  flapping 
cycle.  Details  will  appear  in  [4],  but  are  here  included  because  they  are  not  yet  widely  available. 


3.1  Wing  Trajectories 

An  initial  study  was  directed  at  narrowing  the  pool  of  candidate  kinematic  representations  through 
simulation.  Kinematic  schemes  were  sought  that  were  capable  of  representing  the  most  general  flapping 
patterns,  including  those  with  a  biological  basis,  with  the  fewest  kinematic  parameters.  Candidate 
representations  included  standard  sinusoids,  polar  and  Cartesian  Fourier  series,  piecewise  linear 
segments,  piecewise-defined  sinusoids,  NURBS,  and  scaled  sinusoids. 

Based  on  the  initial  simulation  studies,  the  polar  Fourier  series  and  the  scaled  sinusoids  from 
Berman  and  Wang  (2007)  were  selected  for  further  study.  Using  the  Fourier  series,  the  wing  angles  are 
defined  by 
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where  9S  is  the  sweep  angle,  9V  is  the  pitch  angle,  and  9d  is  the  deviation  angle  of  the  wing;  An  are  the 
kinematic  coefficients;  <pn  are  the  phase  delay  parameters;  and  t'  is  the  nondimensionalized  time  given 
by 

t’  =  t/T  =  ft 

where  /  is  the  flapping  frequency,  T  is  the  period,  and  t  is  the  time.  The  11  parameters  of  the  polar 
Fourier  method  to  be  found  through  hardware-in-the-loop  optimization  are  shown  in  Table  1. 


Table  1.  Polar  Fourier  Kinematic  Parameters 


Input  Parameters 

Symbol 

Stroke  Plane  Angle 

'P 

Sweep: 

Amplitude  of  First  Frequency 

A1.S 

Pitch: 

Constant  Pitch  Offset 

Ao,p 

Amplitude  of  First  Frequency 

ai,p 

Phase  Delay  of  First  Frequency 

01  ,P 

Amplitude  of  Second  Frequency 

A2  ,P 

Phase  Delay  of  Second  Frequency 

02,  P 

Deviation: 

Constant  Pitch  Offset 

ao,d 

Amplitude  of  First  Frequency 

Al  ,D 

Amplitude  of  Second  Frequency 

A2,D 

Phase  Delay  of  Second  Frequency 

02,  D 

A  novel  kinematic  method  was  developed,  based  on  the  pitch  and  sweep  representation  employed 
by  Berman  and  Wang  (2007),  with  two  key  features: 


1.  A  non-sinusoidol  form  of  the  periodic  function.  The  work  resulted  in  the  unification  of  the 
Berman  and  Wang  equations  through  coordinate  transformations  to  enable  continuous 
variation  of  the  shape  from  triangular  to  square  with  the  variation  of  a  single  parameter. 

2.  A  novel  method  of  modifying  the  down-stroke-to-up-stroke  ratio  with  a  single  parameter.  This  is 
done  without  sacrificing  the  continuity  of  the  function  or  its  first  three  derivatives.  This  enables 
the  robotic  mechanism  to  execute  trajectories  with  unequal  half-stroke  durations. 


In  the  Modified  Berman  and  Wang  method  the  wing  angle  trajectories  are  calculated  using 
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Here,  0Sin(t  ,  <Pn,m>  Kn,m)n  an^  0cos(^  >  tpnw  K n,m)n  are  given  by 
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The  relationship  between  anrn  and  Knm  is  given  by 


n  -  1Qa/ioo 


-oo  <  a  <  0,  (Triangular) 


K  =  <  0 


a  =  0,  (Sinusoidal) 


vioa/100  —  1  0  <  a  <  oo,  (Square). 


Additionally,  the  function  r(t')  is  a  "pseudo-time"  function  that  enables  variation  of  the  down- 
stroke-to-up-stroke  ratio.  The  subscript  m  designates  either  sweep  (5),  pitch  (p),  or  deviation  from  the 
stroke  plane  (d).  The  17  parameters  of  the  Modified  Berman  and  Wang  method  to  be  found  through 
hardware-in-the-loop  optimization  are  shown  in  Table  2. 


Table  1:  Inputs  for  Modified  Berman  and  Wang  Kinematics 


Input  Parameters 

Symbol 

Stroke  Plane  Angle 

V 

Down-stroke-to-Up-stroke  Ratio 

D  UratlQ 

Sweep: 

Amplitude  of  First  Frequency 

Ai,s 

Shape  of  First  Frequency 

«i.s 

Pitch: 

Constant  Pitch  Offset 

Ao.p 

Amplitude  of  First  Frequency 

Ai  ,p 

Phase  Delay  of  First  Frequency 

0i,  p 

Shape  of  First  Frequency 

ai  ,P 

Amplitude  of  Second  Frequency 

^2,P 

Phase  Delay  of  Second  Frequency 

02.  P 

Shape  of  Second  Frequency 

a2,P 

Deviation: 

Constant  Pitch  Offset 

Ao,d 

Amplitude  of  First  Frequency 

Ai  ,D 

Shape  of  First  Frequency 

al  ,D 

Amplitude  of  Second  Frequency 

^2  ,D 

Phase  Delay  of  Second  Frequency 

02,  D 

Shape  of  Second  Frequency 

a2,D 

3.2  Experimental  Approach 

The  two  kinematic  representations  described  above  were  used  in  hardware-in-the-loop  experiments 
to  determine  the  optimal  parameters  (Table  1  and  Table  2),  and  therefore  the  optimal  trajectory,  to 
maximize  lift.  Note  that  other  cost  functions  could  be  explored  that  will  optimize  thrust,  lift/thrust 
combinations,  efficiency,  or  other  flight  parameters.  The  optimization  method  is  similar  to  a  Box- 
Behnken/Response  Surface  method,  but  with  modifications  that  include  intelligent  step  size  selection, 
exclusion  of  insignificant  parameters  from  the  kinematic  models,  and  termination  criteria. 

The  same  hardware  was  used,  with  the  exception  that  a  new  shear  force  sensor  was  incorporated 
into  the  wing,  allowing  a  more  accurate  measurement  of  wing  forces.  The  wing  angles  were  defined 


using  the  sweep,  pitch,  and  deviation  angles  described  above,  and  converted  into  joint  angles  using  a 
series  of  coordinate  transformations.  Tracking  of  each  joint  was  achieved  with  decoupled  PID  loops 
implemented  in  National  Instruments  real-time  hardware. 

3.3  Optimization  Results 

Initial  conditions  were  assigned  for  each  of  the  kinematic  representations  and  the  modified  Box- 
Behnken  optimization  routine  was  used  to  cover  the  parameter  space  and  determine  the  kinematic 
parameters  that  optimize  lift  on  the  robotic  flapping  mechanism.  The  evolution  of  the  optimization  for 
both  kinematic  representations  is  shown  in  Fig.  3.  Figure  3a  shows  a  plot  of  lift  vs.  thrust  for  each 
iteration  of  the  optimization  process.  As  can  be  seen  from  the  figure,  as  the  hardware-in-the-loop 
optimization  proceeded,  trajectories  were  found  that  resulted  in  improved  lift  for  both  representations. 
Figure  3b  shows  the  evolution  of  the  fitness  (cost)  function  as  the  optimization  proceeds,  illustrating 
that  the  method  results  in  decreasing  values  of  the  fitness  function,  implying  that  the  each  iteration 
improves  over  the  previous. 
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Figure  3.  Evolution  of  lift  and  thrust  (left)  and  fitness  function  (right)  for  successive  optimization 

iterations. 


The  evolution  of  the  kinematic  parameters  is  shown  in  Table  3  (polar  Fourier  series)  and  Table  4 
(Modified  Berman  and  Wang).  The  tables  are  color  coded  to  highlight  those  regions  where  significant 
change  is  going  on.  For  the  polar  Fourier  series  (Table  3),  the  most  significant  changes  are  observed  for 
the  amplitude  of  the  sweep  angle  (i4ls)  and  the  amplitude  of  the  second  deviation  angle  harmonic 
(^2  ,d)‘ 


Table  3.  Center  Point  Changes,  Polar  Fourier  Series 


Aid 

Aqj, 

<Pij> 

A2j> 

<p2J> 

Aoj) 

A\j> 

Aid 

<t>2D 

Iteration  2 

-101 

15 

0.80 

-9.24 

-3.95 

0.50 

0 

2.96 

-1.50 

-9.32 

0 

Iteration  3 

0 

8.56 

0 

0 

0 

0 

0 

0 

-7.98 

-8.42 

15 

Iteration  4 

0 

6.27 

0 

0 

0 

0 

0 

0 

3.82 

-4.00 

0 

Iteration  5 

0 

4.64 

0 

0 

0 

0 

0 

0 

4.38 

-4.72 

-15 

Iteration  6 

0.56 

3.40 

-3.58 

-4.52 

7.39 

5.46 

0 

3.34 

-0.70 

-3.67 

5.43 

Optimum 

0.23 

0.96 

-1.37 

-0.99 

0.45 

4.48 

9.70 

0.39 

0.57 

-0.64 

228 

Table  2.  Center  Point  Changes,  Modified  Berman  anc 


Wang 


V 

W ratio 

Ais 

<*IS 

Au> 

<Pi  J> 

J> 

AZJ> 

<PlJ> 

a2j> 

Aoa 

aU> 

A  2  J> 

<Pzj) 

aZD 

Iteration  2 

0 

20 

12.57 

0 

0 

-7.39 

0 

-30 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Iteration  3 

-0.71 

-7.73 

9.32 

8.19 

-1.40 

-9.17 

-0.62 

B.39 

1.15 

0 

0 

3.70 

-0.88 

0 

-6.74 

0 

0 

Iteration  4 

0 

0 

13.65 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-15 

0 

0 

Iteration  5 

0 

0 

6.41 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-9.46 

0 

0 

Optimum 

0.79 

-111 

0.20 

-0.65 

-0.99 

IB 

1.73 

10.66 

1.54 

0 

0 

-4.12 

-0.38 

0 

-4.70 

0.59 

0 

For  the  Modified  Berman  and  Wang  method  (Table  4),  it  becomes  evident  here  that  certain 
parameters  begin  the  optimization  routine  with  large  step  sizes  which  are  reduced  for  later  iterations. 
Again,  the  primary  example  of  this  is  the  sweep  amplitude  (711S).  These  parameters  tend  to  dominate 
the  effect  on  the  fitness  function  early  on,  but  that  influence  is  reduced  as  the  wing  motion  approaches 
the  limits  of  the  joint  angles. 

Other  observations  can  be  made  from  this  data  as  well.  For  instance,  there  are  variables  unique  to 
the  Modified  Berman  and  Wang  trajectory  generator  which  are  not  used  at  all  during  the  test  (a2P, 
a1D,  oc2r d).  Of  the  parameters  unique  to  Modified  Berman  and  Wang,  three  moved  through  the  design 
space.  The  first  of  these,  the  down-stroke-to-up-stroke  ratio  modifier  (DUratio),  was  raised  by  20%,  the 
full  value  allowed  during  the  first  iteration.  In  the  second  iteration  it  returned  back  toward  its  origin, 
ending  at  12.27%.  During  the  final  iteration,  it  adjusted  downward  slightly  to  finish  at  11.15%.  This 
means  that  for  the  initial  portion  of  the  optimization  routine,  benefit  was  derived  by  lengthening  the 
duration  of  the  down-stroke,  while  later  in  the  process,  benefit  was  had  by  equalizing  the  durations  of 
the  two  half  strokes.  The  other  active  parameters  were  the  shape  of  the  sweep  and  pitch  fundamental 
frequencies  (als  and  a1P). 

Parameters  not  unique  to  Modified  Berman  and  Wang  were  also  fairly  inactive  throughout  the 
process.  Table  2  shows  that  these  included  the  pitch  offset  (A0P),  the  phase  shift  of  the  first  pitch 
harmonic  (0 1P),  the  second  amplitude  in  pitch  (A2P),  and  the  second  phase  shift  in  pitch  (02;P)  among 
others.  Additionally,  the  stroke  plane  angle  W  stayed  fairly  constant  throughout  the  process. 

4.  EFFICIENT  HARDWARE-IN-THE-LOOP  OPTIMIZATION 

Hardware-in-the-loop  (HIL)  modeling  is  the  use  of  physical  hardware  to  represent  a  system  (as 
opposed  to  an  analytical  or  computational  model).  One  advantage  of  HIL  modeling  is  that  it  can  avoid 
problems  caused  by  inaccurate  models.  Complicated  phenomena  may  be  difficult  to  model  accurately, 
but  they  can  act  on  a  hardware  model  similarly  to  how  they  would  act  on  a  real  system.  Another  benefit 
of  HIL  modeling  is  manifest  when  accurate  models  cannot  be  obtained  because  of  time  or  other 
constraints.  In  cases  like  these,  a  system  can  be  modeled  by  a  HIL  model,  which  is  treated  as  a  black  box 
that  returns  performance  outputs  given  inputs. 

Some  researchers  have  used  HIL  in  the  study  of  flapping  flight.  However,  HIL  has  some  challenges 
when  used  in  optimization.  Sometimes,  HIL  models  are  too  expensive  to  call  exhaustively.  The  HIL 
system  used  previously  [1,2]  has  a  large,  complex  design  space  and  sampling  the  whole  space  with  a 
central  composite  design  of  experiments  could  require  over  32,000  calls.  While  this  test  would  take  a 
long  time  to  complete,  it  would  also  result  in  system  wear.  In  cases  like  this,  mechanical  wear  can 
contaminate  the  results  of  the  test.  This  makes  the  hardware  too  expensive  to  use  in  a  normal 
optimization  algorithm.  Instead,  a  less  expensive  model  is  needed  for  optimization.  In  this  research,  the 
use  of  a  physics-based  model  as  a  surrogate  was  studied  as  a  means  for  reducing  the  number  of 
hardware  calls  while  maintaining  the  benefits  of  HIL  modeling.  The  below  is  a  summary  of  the  methods 
and  results  presented  in  [5,6]. 


4.1  Approach 

An  approach  to  variable  fidelity,  hardware-in-the-loop  optimization  that  minimizes  calls  to  the 
hardware  is  here  presented.  The  basis  of  this  approach  is  an  intelligent  DOE  (IDOE)  that  is  informed  by 
physics.  As  such,  the  intent  of  this  research  was  to  create  an  IDOE.  An  IDOE  chooses  points  to  test  that 
are  predicted  to  have  better  objectives  than  the  previous  points.  This  requires  something  that  can 
provide  accurate  information  about  the  hardware  without  actually  calling  the  hardware.  A  physics-based 
surrogate  model  can  fill  this  role.  The  physics-based  model  predicts  optimal  solutions  and  gives 
information  about  how  different  design  variables  affect  the  performance  of  the  hardware.  This 
information  can  be  used  to  choose  the  next  point  or  points  for  the  DOE  to  test. 

4.2  Intelligent  DOE 

In  order  to  take  advantage  of  the  similarities  between  a  physics-based  model  and  a  hardware 
model,  some  mapping  between  the  expensive  hardware  model  and  the  inexpensive  physics  model  must 
be  made.  This  mapping  allows  the  physics  model  to  provide  useful  information  about  the  hardware 
model  by  relating  points  within  the  physics  model  space  to  corresponding  ones  in  the  hardware  model 
space.  In  order  to  create  this  mapping,  some  assumptions  must  be  made  about  the  differences  between 
the  models.  In  this  research,  the  assumption  was  that  a  point  in  the  hardware  model  space,  xh,  can  be 
represented  as  a  point  in  the  physics  model  with  some  shift,  y  Also,  the  objective  values  may  differ  by 
some  scaling  factor  a.  The  two  models,  then,  can  be  related  by  the  expression  H(xh)  =  aP(xp  -  fj,  where 
xp  is  a  point  within  the  physics  model  space  and  xh  is  a  point  within  the  hardware  model  design  space.  H 
represents  the  hardware  model  and  P  represents  the  physics  model.  Since  xp  is  an  element  within  the 
design  space  and  y  is  defined  within  the  same  space,  both  are  vectors  of  the  same  number  of 
dimensions  as  the  physics  model  has.  There  are  other  potential  differences  between  the  models,  such  as 
stretching;  however,  for  many  situations,  the  differences  discussed  here  are  sufficient  to  map  the 
models  sufficiently  well  for  optimization  purposes. 

4.3  Algorithm 

Here,  the  steps  in  the  intelligent  DOE  are  outlined.  The  optimization  problem  can  be  stated  as 

min  H(x) 


subject  to 


linear  constraints  <  0 
nonlinear  constraints  <  0 


where  H(x)  is  the  hardware  model. 

The  next  step  to  finding  the  hardware  model  solution,  jca*,  is  to  find  values  for  a  and  y  for  which 
H(x*)  =  aP(xp /  -  y),  where  xp*  =  xp/  -  y,  and  xpf  is  the  initial  optimal  solution  of  the  physics  model. 
When  these  values  have  been  found,  Xh  =xp  . 


; — Model  PH* 


Begin 


Find  fa 


Figure  4.  Flowchart  of  the  intelligent  DOE 


Figure  4  shows  a  flowchart  of  the  entire  method.  First,  the  user  finds  the  optimal  solution,  xp* ,  to 
the  low  cost  physics  model.  This  involves  setting  bounds  on  the  space,  choosing  a  starting  point,  and 
using  a  known  optimization  technique.  The  optimization  problem  is  expressed  by 


min  P(xp) 

Xp 


subject  to 


linear  constraints  <  0 
nonlinear  constraints  <  0 


where  P(pc )  is  the  physics  model. 

In  this  research,  gradient-based  algorithms  were  used,  however,  genetic  or  other  algorithms  could 
be  effective,  depending  on  the  nature  of  the  model.  is  the  objective  value  at  xp/.  Step  2  sets  xp*  = 
Xpj  .  xp>i  is  the  initial  solution  to  the  physics  model  and  does  not  change  throughout  the  algorithm,  while 
xp  does  change  as  the  physics  model  is  shifted.  Step  3  calls  the  hardware  model,  H,  at  the  physics 
model  solution,  xp\  and  sets  the  result  as  jui,.  Step  4  finds  the  gradient  of  the  hardware  model  at  xp  . 
Steps  5  through  10  are  repeated  until  the  algorithm  converges.  Once  VH(xp)  has  been  found,  step  5 
finds  the  point,  xP,„,  where  the  physics  model  gradient  matches  the  gradient  measured  in  the  hardware. 
This  step  is  an  optimization  problem  that  can  be  stated  as 

minL(x) 


subject  to 


linear  constraints  <  0 
nonlinear  constraints  <  0. 


Where 


L=WVH(x;)-VP(xp.„)\\ 


This  is  done  with  known  optimization  techniques.  Here,  gradient-based  and  genetic  algorithms  were 
used  as  needed  by  the  problem.  Step  6  finds  the  offset  between  the  hardware  and  physics  models,  /  by 
subtracting  xp>m  from  xp *.  This  difference  is  added  to  the  existing  offset,  which  is  initially  zero,  and  the 
physics  model  is  shifted  by  this  amount.  Next,  step  7  defines  the  new  xp  by  shifting  it  by  the  offset  so 
that  it  does  not  move  relative  to  the  physics  model.  /4  and  VH are  found  at  this  new  xp  in  step  8.  Step  9 
updates  a  by  dividing  the  best  value  returned  by  the  hardware,  nh,b,  by  nP;i,  the  value  of  the  physics 
model  at  xp/.  This  scale  must  be  accurate  in  order  for  the  physics  gradients  to  be  correct.  Finally,  step 
10  checks  for  convergence.  This  can  be  done  using  any  appropriate  criteria.  A  convergence  criteria  for 
variable  fidelity  optimization  is  presented  in  the  literature,  however,  this  uses  a  traditional  response 
surface,  which  allows  for  many  calls.  This  method  can  require  a  prohibitive  number  of  function  calls 
when  calling  the  expensive  physics  model.  Two  alternative  convergence  criteria  requiring  no  additional 
hardware  calls  are  presented  here.  One  criterion  is  whether  the  gradient  was  zero.  As  is  the  case  with 
some  gradient-based  algorithms,  this  can  terminate  the  algorithm  at  local  maxima,  minima,  or  saddle 
points,  so  care  should  be  taken  to  ensure  that  solutions  reached  with  these  criteria  are  optimal.  In  many 
cases,  solutions  to  engineering  problems  are  constrained,  which  means,  among  other  things,  that  the 
gradients  are  non-zero.  For  cases  like  this,  convergence  can  be  triggered  when  VP  (xp*)  =  VH  (xh*).  This 
criterion,  in  effect,  is  triggered  when  the  hardware  point  that  corresponds  to  the  physics  model  solution, 
xp,  has  been  found.  This  assumes  that  the  scale,  a,  and  the  offset,  /have  been  found  correctly.  Other 
convergence  criteria  could  be  used,  depending  on  the  nature  of  the  problem  being  solved.  Using  these 
steps,  the  physics  model  shifts  and  scales  until  it  matches  the  hardware  model  at  the  solution  and  xp  = 
xh*.  As  can  be  seen  in  the  flowchart,  each  time  steps  5-10  are  iterated,  the  hardware  gradient  must  be 
found.  Thus,  the  number  of  hardware  calls  needed  for  the  IDOE  to  converge  is  directly  related  to  the 
number  of  iterations  of  the  algorithm  and  it  is  desirable  to  find  the  solution  with  as  few  iterations  as 
possible. 

Figure  5  shows  plots  of  hypothetical  hardware  and  physics  models  after  each  iteration  of  a  one¬ 
dimensional  test.  The  objective  of  this  problem  is  to  minimize  the  value  returned  by  the  hardware.  Here, 
the  hardware  model  has  been  scaled  by  a  factor  of  a=  1.3  and  shifted  by  y=  20°.  Figure  5a  shows  the 
models  in  their  initial  state.  Notice  that  it  has  larger  negative  value  than  which  is  caused  by  the  fact 
that  the  hardware  model  is  scaled  larger  than  the  physics  one.  In  iteration  2,  a  and  /have  been  updated 
so  that  the  physics  model  is  nearly  equal  to  the  hardware  model.  By  the  third  iteration,  shown  in  Fig.  5c, 
the  models  are  equal  and  the  solution  to  the  hardware  model,  xh*  equals  the  known  solution  to  the 
physics  model,  xp. 
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Figure  5.  Progression  of  the  IDOE  through  a  sample  optimization  routine. 

4.4  Validation  and  Verification 

After  development,  the  intelligent  DOE  (IDOE)  was  verified  on  several  problems  of  varying 
complexity.  This  allowed  for  testing  of  a  variety  of  scenarios.  Each  example  problem  was  optimized  by 
the  IDOE  and  two  DOE-based  VFO  routines  for  comparison.  The  results  are  provided  in  [5]. 

5.  SYNTHETIC  APERTURE  PARTICLE  IMAGE  VELOCIMETRY 

Three-dimensional  velocity  measurements  of  the  flow  fields  around  a  tethered  mechanical  flapper 
were  obtained  using  synthetic  aperture  PIV  (SAPIV).  The  results  focused  mainly  on  the  leading  and 
trailing  edge  vortices  (LEV,  TEV)  through  an  entire  flapping  cycle.  The  results  were  compared  to  velocity 
measurements  taken  using  2D  PIV.  Additionally,  force  measurements  of  the  lift  and  thrust  generated  by 
the  mechanical  flapper  were  taken.  The  visual  hull  of  the  mechanical  flapper  was  reconstructed.  The 
methods  and  results  are  being  prepared  for  a  manuscript  [7]  (anticipated  for  submission  Fall  2013),  but 
are  included  here  since  they  are  not  yet  widely  available. 

Measurements  were  performed  on  a  commercially  available  MAV,  powered  by  a  small  electric 
motor  that  drives  a  pivoting  mechanism  at  the  front  edge  of  the  wings,  with  the  back  edge  of  the  wings 
pinned  (Fig.  6).  A  strut  forms  the  leading  edge  of  the  wings  with  the  rest  of  the  wings  made  of  a  thin, 
plastic  membrane.  Steering  is  accomplished  by  moving  the  pinned  back  edges  of  the  wings  up  or  down, 


depending  on  the  direction  of  the  turn.  The  mechanism  is  controlled  via  wireless  remote.  Stability  is 
achieved  through  a  fixed  tail  section  connected  to  the  back  of  the  body.  The  wing  beat  frequency  was 
15.9  Hz,  resulting  in  a  Reynolds  number  (calculated  using  the  induced  forward  velocity  and  mean  chord 
length)  of  9600  and  a  Strouhal  number  of  0.79. 

The  experiments  were  performed  in  a  61  x  41  x  30  cm3  custom  acrylic  observation  tank,  providing 
complete  optical  access  (Fig.  6).  The  MAV  was  tethered  to  a  steel  plate  in  the  center  of  the  tank.  The 
synthetic  aperture  setup  consisted  of  eight  Photron  SA3  high-speed  video  cameras.  The  cameras  were 
equipped  with  50  mm  lenses  and  arranged  in  an  array  such  that  every  camera  was  looking  at  the  same 
region  of  interest.  The  camera  array  was  positioned  parallel  to  the  long  axis  of  the  observation  tank.  A 
Darwin-Duo  laser  system  (Quantronix,  Nd:YLF,  1000  Hz)  was  used  to  provide  PIV  illumination.  This  laser 
uses  two  laser  pulses  at  500  Hz  each,  corresponding  to  PIV  at  500  Hz.  The  laser  was  positioned 
perpendicular  to  the  cameras.  A  beam  expander  from  Edmunds  Optics  enlarged  the  beam  to  a  cross 
sectional  area  of  55.8  cm2.  The  air  inside  the  tank  was  seeded  with  hollow  polymer  microspheres  with  a 
mean  diameter  of  40  |um  and  a  density  of  25  kg/m3  (Expancel,  461  DET  40  d25). 


Figure  6.  Line  drawings  of  the  MAV  (left)  and  experimental  setup  (right).  Eight  high-speed  cameras  were 
used  to  perform  SAPIV.  A  ninth  camera,  above  the  tank,  could  be  used  for  wing  tracking. 


Three-dimensional,  time-resolved  fluid  velocities  were  obtained  by  using  3D  SAPIV.  SAPIV  is  based 
on  the  principles  of  light  field  imaging.  The  essential  idea  is  that  an  array  of  cameras  can  image  an  object 
from  multiple  viewpoints.  Since  each  viewpoint  sees  the  object  from  a  different  angle,  the  images  from 
each  camera  can  be  digitally  refocused  to  create  a  set  of  images  that  are  in  focus  at  different  planes 
throughout  the  region  of  interest,  known  as  a  focal  stack.  The  use  of  multiple  viewpoints  has  several 
advantages,  such  as  allowing  for  high-density  seeding,  because  of  the  ability  to  see  around  partial 
occlusions.  Particles  that  are  in  focus  on  a  given  plane  will  have  a  higher  intensity  than  particles  that  are 
not  in  focus.  The  out-of-focus  particles  can  be  eliminated  from  the  image  planes  using  thresholding.  The 
digitally  refocused  and  thresholded  image  pairs  can  then  be  passed  to  a  standard  3D  PIV  system  such  as 
matPIV  for  cross  correlation  (Belden  et  al,  2010).  The  cameras  were  calibrated  using  a  process  based  on 
the  Multi-camera  Self  Calibration  Method  developed  by  Svoboda  et  al.  (2005).  The  preprocessed  images 
from  each  of  the  eight  cameras,  and  the  camera  projection  matrix,  were  passed  to  the  SAPIV  software.  A 


map-shift-average  algorithm  was  used  to  digitally  refocus  the  images  onto  synthetic  focal  planes  to 
create  a  focal  stack.  SAPIV  code  is  available  at  www.saimaging.org. 

The  lift  and  thrust  forces  generated  by  the  flapping  wing  of  the  MAV  were  measured  with  a  single 
axis  load  cell.  Force  measurements  were  taken  for  10  different  runs,  with  55  flapping  periods  in  each  run 
for  statistical  analysis.  The  MAV  was  mounted  in  two  configurations  -  horizontally  for  lift  measurements, 
and  vertically  for  thrust  measurements.  A  single  SA3  camera  was  synchronized  with  the  load  cell  to 
document  the  position  of  the  MAV  wings  as  force  data  points  were  acquired.  The  load  cell  sampled  at 
4000  Hz  and  the  camera  sampled  at  2000  Hz. 

Figure  7  shows  the  wing  of  the  MAV  through  a  partial  (last  75%)  down-stroke  in  the  left  two 
columns  and  a  partial  (first  75%)  up-stroke  in  the  right  two  columns,  with  vorticity  plotted  as  iso¬ 
surfaces.  The  presence  of  an  LEV  and  a  TEV  are  recognized  on  the  leading  edge  (left  side)  and  trailing 
edge  (right  side)  of  the  wing,  respectively.  The  wing  was  recreated  in  the  plots  using  the  visual  hull 
method.  The  recreated  wing  is  not  visible  in  the  plots  when  the  actual  wing  drops  below  the  laser 
volume  and  at  times  when  the  amount  of  particles  gives  the  visual  hull  algorithm  difficulty.  However, 
LEVs  and  TEVs  are  still  visible  in  their  respective  positions  in  the  plots. 

The  LEV  and  TEV  remain  attached  to  the  wing  through  the  entire  down-stroke  and  up-stroke.  The 
vortices  are  shed  into  the  wake  during  the  transitions  between  the  up  and  down  strokes.  This  is  evident 
at  time  28  ms  when  the  wing  is  transitioning  from  down-stroke  to  up-stroke.  The  LEV  from  the  down- 
stroke  on  the  bottom  side  of  the  wing  (center  structure)  is  seen  leaving  the  leading  edge  and  moving 
into  the  wake  as  a  new  LEV  (left  structure)  is  seen  forming  on  the  leading  edge  on  the  top  side  of  the 
wing.  During  the  down-stroke,  the  LEV  and  TEV  form  and  remain  attached  to  the  top  of  the  wing,  and 
likewise  during  the  up-stroke,  the  LEV  and  TEV  form  and  remain  attached  to  the  bottom  of  the  wing.  A 
top  view  of  the  wing  shows  the  varying  shape  of  the  LEVs  and  TEVs  along  the  span  of  the  wing  (Fig.  10). 
A  look  at  the  cross  section  shows  that,  at  mid-span  of  the  wings,  the  LEV  has  a  minimum  diameter.  The 
LEV  shows  greater  size  at  the  wing  tip. 

Figure  9  shows  vorticity  slices  from  the  3D  velocity  fields  at  65%  of  the  half  span,  as  well  as  the 
vorticity  plots  from  the  2D  PIV  at  60%  half-span,  at  14  time  steps  from  mid-down-stroke  through  an 
entire  flapping  period  (8  ms  to  60  ms).  The  comparison  is  promising,  showing  agreement  in  the 
formation  and  behavior  of  the  LEV  and  TEV  on  the  MAV  wing.  In  both  cases,  the  LEV  is  seen  to  remain 
attached  to  the  leading  edge  of  the  wing  throughout  the  down-stroke.  The  growth  of  the  LEV  can  also 
be  seen  through  the  four  time  steps  shown  of  the  down-stroke.  In  the  last  frame  the  wing  enters  the 
transition  between  down-stroke  and  up-stroke.  The  shedding  of  the  TEV  is  visible  and  the  beginning  of 
the  shedding  of  the  LEV  is  also  visible. 

Differences  are  seen  in  the  flow  beneath  the  wing.  In  the  2D  PIV  plots,  vortices  are  seen  under  the 
wing  that  are  not  visible  in  the  SAPIV  slices.  This  is  due  to  the  2D  PIV  plots  being  shown  after  a  previous 
flapping  period  while  the  SAPIV  slices  are  from  the  first  flapping  period.  The  vortices  in  the  2D  PIV  plots 
are  the  shed  LEV  and  TEV  from  the  up-stroke  previous  to  the  down-stroke  shown  in  the  plots.  It  appears 
that  the  vortical  structures  from  the  previous  stroke  are  partially  recaptured  in  the  current  stroke 
instead  of  being  entirely  shed  into  the  wake. 

The  vortical  structures  in  the  3D  PIV  plots  are  more  apparent  in  the  down-stroke  than  the  up-stroke. 
This  is  due  to  the  turbulence  introduced  inside  the  acrylic  observation  tank  as  the  MAV  wings  move  the 
air  around.  The  up-stroke  comes  after  the  initial  down-stroke  in  the  experiment;  thus,  the  LEV  and  TEV 
in  the  up-stroke  are  less  distinguishable  than  in  the  down-stroke. 


Figure  7.  Iso-surfaces  of  vorticity  about  the  z-axis  plotted  for  12  time  steps  covering  0.75  wing  beats.  A 
down-stroke  persists  from  1  ms  to  21  ms.  50%  of  the  up-stroke  is  shown  from  25  ms  to  45  ms. 
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Figure  8.  A  comparison  of  the  SAPIV  (left  columns)  and  2DPIV  (right  columns)  data  with  the  applying 
period  shown  in  14  time  steps. 
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Figure  9.  Top  view  of  the  wing  at  mid-down-stroke.  The  diameter  of  the  LEV  (left)  and  TEV  (right)  vary 
from  the  wing  mid  span  (top)  to  the  wing  tip  (bottom). 


The  visual  hull  of  the  wing  represents  the  volume  in  which  the  wing  is  contained,  but  requires  large 
connected  components  for  the  algorithm  to  process.  The  reconstruction  of  the  wing  could  be  improved 
by  less  particle  seeding  or  by  adding  more  light  sources  so  that  the  intensity  of  the  wing  is  the  same  in 
all  camera  views  to  assist  in  connected  component  identification.  In  the  results  presented  in  Fig.  7,  the 
MAV  starts  with  its  wing  at  an  angle  of  about  45°  to  the  z-axis  but  as  the  flight  progresses  the  wing 
reaches  an  angle  of  about  -45°  with  respect  to  the  z-axis.  By  using  multiple  cameras,  particles  that  are 
obscured  in  one  camera  view  are  visible  in  other  camera  views,  thus  providing  whole-field 
measurements  at  each  wing  position.  The  2D  PIV  vorticity  plots  show  greater  detail  in  the  LEV  and  TEV. 
Despite  the  lower  resolution  in  the  3D  SAPIV  data,  it  is  possible  to  see  similarities  in  the  structure  and 
behavior  of  the  vortices.  Greater  resolution  for  the  SAPIV  results  can  be  attained  by  reducing  the  size  of 
the  volume  of  interest.  This  will  make  it  more  difficult  to  visualize  the  entire  flapping  period  in  one 
experimental  run,  but  would  provide  greater  detail  for  selected  phases  of  the  flapping  period. 

The  wing  position  was  determined  by  tracking  the  wing  tip  of  the  MAV  and  plotting  the  pixel 
locations.  It  was  observed  that  the  MAV  has  an  interesting  wing  tip  trajectory  in  the  x-y  plane  (Fig.  10a). 
The  wing  is  driven  by  an  electric  motor  through  a  four-bar  linkage,  providing  an  up  and  down  motion. 
However,  the  wing  tip  follows  a  skewed  figure-eight  trajectory,  with  the  figure-eight  skewed  in  the 
forward  direction. 

A  further  investigation  of  the  wing  tip  trajectory  was  performed  by  including  a  second  SA3  camera  in 
the  acquisition  loop  and  performing  a  direct  linear  transform  (DLT)  analysis  of  the  wing  tip  position.  The 
trajectory  of  the  wing  tip  in  the  x-y,  x-z,  and  y-z  planes  can  be  seen  in  Fig.  10.  The  skewed  figure-eight 
trajectory  is  seen  in  the  x-y  plane  (Fig.  10a),  and  the  wing  trajectory  has  a  curvature  with  a  radius  equal 
to  the  span-wise  length  of  the  wing  (Fig.  10c).  The  trajectory  in  the  x-z  plane  (Fig.  10b)  appears  to 
deviate  from  the  symmetric  trajectory  pattern  as  in  the  other  planes;  however,  it  can  be  seen  in  Fig.  10c 
that  the  wing  is  not  flapping  symmetrically  with  the  x,  y,  and  z  axes  chosen.  However,  the  same 
symmetry  is  exhibited  along  an  axis  skewed  from  the  chosen  axes. 

Figure  11  shows  the  mean  lift  and  thrust  forces  generated  from  55  flapping  periods.  The  force  data 
was  processed  using  a  3rd-order,  120  Hz  low-pass  Butterworth  filter  to  remove  a  150  Hz  signal  present  in 
both  the  lift  and  thrust  force  measurements  generated  by  the  motor  and  gearing  of  the  four-bar  linkage 
driving  the  wings.  The  plot  is  split  down  the  center  to  distinguish  the  up-stroke  from  the  down-stroke  of 
the  flapping  period. 

In  the  up-stroke,  the  lift  force  is  negative,  with  a  minimum  near  the  end  of  the  up-stroke  occurring 
at  -0.58  N.  The  primary  generation  of  the  negative  lift  force  is  the  downward  force  generated  by  the 
upward  movement  of  the  wing.  There  is  also  negative  lift  generated  by  the  LEV  witnessed  on  the 
underside  of  the  wing  during  the  up-stroke.  During  the  down-stroke,  the  forces  are  reversed;  an  upward 
force  is  generated  by  the  downward  movement  of  the  wing,  and  positive  lift  is  generated.  The  maximum 
lift  force  occurs  during  the  second  half  of  the  down-stroke  at  0.6405  N.  The  lift  is  generated  by  the  LEV 
present  on  top  of  the  wing  during  the  down-stroke. 
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Figure  10.  3-dimensional  trajectory  of  the  wing  tip  as  determined  by  performing  a  DLT  analysis.  The 
trajectory  is  plotted  on  the  XY  (a),  XZ  (b),  and  Y  Z  (c)  planes. 


Average  Fence  Measurements 


Figure  11.  Lift  and  thrust  measurements  taken  using  a  single-axis  load  cell.  The  plots  represent  the 
forces  seen  over  one  flapping  period,  averaged  over  55  flapping  periods.  Error  bars  are  plotted  for  the 
lift  and  thrust  measurements. 

The  weight  of  the  MAV  was  0.1275  N  (calculated  using  the  load  cell).  By  taking  the  average  force 
over  the  entire  flapping  period,  the  net  lift  and  net  thrust  (less  the  weight  of  the  MAV)  were  determined 
to  be  0.014  N  and  0.084  N,  respectively,  with  95%  confidence.  This  translates  to  a  net  lift  force  11.0% 
greater  than  the  weight  of  the  MAV  and  a  net  thrust  force  65.9%  greater  than  the  weight  of  the  MAV. 

The  vortical  structures  on  the  MAV  were  observed  using  both  SAPIV  and  2DPIV.  Changing  cross 
sections  of  the  LEV  and  TEV  show  evidence  of  spanwise  flow  in  the  direction  of  the  wing  tip.  The  3D 
vorticity  plots  and  2D  slice  vorticity  plots  show  the  three-dimensional  nature  of  the  LEV  and  TEV.  And 
finally,  the  2DPIV  plots  and  the  2D  slices  from  the  3D  data  show  general  agreement  in  the  structure  and 
behavior  of  the  flow  around  the  flapping  wing.  Further  work  can  be  done  to  study  the  expected  span 
wise  flow  using  smaller  search  volumes  and  larger  experimental  volumes  to  reduce  air  circulation  within 
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the  acrylic  box. 

Multiple  viewpoints  used  in  SAPIV  allow  for  the  ability  to  see  around  partial  occlusions,  including 
high  seeding  densities.  SAPIV  is  able  to  reconstruct  3D  whole-field  velocity  fields  in  large  spatial  regions 
of  interest  and  the  whole-field  nature  of  the  velocity  fields  allows  for  analysis  of  flow  characteristics  not 
rectified  with  the  camera  axes.  The  results  demonstrate  that  SAPIV  can  be  used  to  measure  fluid  flow 
velocities  and  reconstruct  the  visual  hull  of  flapping  wings,  with  anticipated  application  in  studying  the 
complex  and  unsteady  nature  of  this  flight  regime. 
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